
#include <math.h>
#include <stdio.h>

#include "Component.h"
#include "StaticDeformationState.h"
#include "FEMModel.h"
#include "TopologicalOperation.h"
#include "Tetrahedron.h"
#include "TetrahedronState.h"
#include "LinearCGState.h"
#include "BigVector.h"
#include "NonlinearCGState.h"
#include "Node.h"
#include "Statistics.h"
#include "IncrementalForceIteration.h"
#include "IncrementalNewtonSolver.h"

#include "Types.h"

StaticDeformationState::StaticDeformationState (FiniteElementModel*fem)
: DeformationState (fem)
{
	iter_state_ = 0;

	bool le =!strcmp (GetStringSetting ("elasticity"), "nonlinear");
	char const * relt = GetStringSetting ("relaxation-type");
	if (!strcmp (relt, "linear"))
	{
		if (le)
			fprintf (stderr, " **** Requested non-linear elasticity and linear CG. Using nonlinear CG instead.;\n");
		else
			iter_state_ = new LinearCGState ();
	}
	else if (!strcmp (relt, "incremental"))
	{
		iter_state_ = new IncrementalForceIteration (this);
	}
	else if (!strcmp (relt, "incremental-newton"))
	{
		iter_state_ = new IncrementalNewtonSolver (this);
	}
	else if (!strcmp (relt, "nonlinear"))
	{
		iter_state_ = new NonlinearCGState ();
	}
	else
	{
		Warning ("No such solver: using nonlinear CG");
	}

	if (!iter_state_)
		iter_state_ = new NonlinearCGState ();
}

void StaticDeformationState::UpdateVectorSizes ()
{
	DeformationState::UpdateVectorSizes ();
	iter_state_->Reinit (this);
}

void StaticDeformationState::UpdateForces ()
{
	DeformationState::UpdateForces ();
	iter_state_->UpdateResidual (this);
}

bool StaticDeformationState::GoodSolution ()const
{
	bool b = iter_state_->GoodSolution (this);

	if (b && iter_state_->iter_count_)
	{
		int it = iter_state_->iter_count_ ;
		printf ("Solution found in %d iterations\n ",  it);
		printf ("Iters / cbrt (nodecount) =%f\n" , it / pow (fem_->NodeCount (), 1.0/3.0));

		iter_state_->iter_count_ = 0;
		if (!quiet_global_b)
			StatisticTick ();	  
	}
	return b;
}


void StaticDeformationState::DoSimulationBody ()
{
	static bool last_time_ok ;
	int warncount  = 100000 / (Dimension () + 1) + 1;  

	if (GoodSolution ())
	{
		if (!last_time_ok)
			printf ("solution reached\n");
		last_time_ok = true;
	}
	else
	{
		last_time_ok = false;
		// if (!iter_state_->degenerate_b_)
		iter_state_->IterationStep (this);
		if (!(global_iteration_count % warncount ))
			StatisticTick ();
	}
}

/*
In a static simulation, every component must be fixed to ascertain a
unique solution.  */
void StaticDeformationState::InitSimulation ()
{
	ComponentMap * cs = FindComponents (fem_);
	FixLooseComponents (fem_, cs);
	delete cs;

	DeformationState::InitSimulation ();
}


void StaticDeformationState::CheckComponents (ComponentMap const*cs)
{
	if (!cs )
	{
		cs = FindComponents (fem_);
	}

	FixLooseComponents (fem_, cs);
}




LinkArray< Array<Real> > StaticDeformationState::VectorVariables ()
{
	LinkArray<Array <Real> > vvs =  DeformationState::VectorVariables ();
	vvs.Concat ( iter_state_->VectorVariables ());
	return  vvs;
}
